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Abstract The functioning of a hving cell is largely de- 
termined by the structure of its regulatory network, 
comprising non-linear interactions between regulatory 
genes. An important factor for the stability and evolv- 
ability of such regulatory systems is neutrality — typ- 
ically a large number of alternative network structures 
give rise to the necessary dynamics. Here we study the 
discretized regulatory dynamics of the yeast cell cycle 
[Li et al., PNAS, 2004] and the set of networks capa- 
ble of reproducing it, which we call functional. Among 
these, the empirical yeast wildtype network is close to 
optimal with respect to sparse wiring. Under point mu- 
tations, which establish or delete single interactions, the 
neutral space of functional networks is fragmented into 
^ 4.7 X 10^ components. One of the smaller ones con- 
tains the wildtype network. On average, functional net- 
works reachable from the wildtype by mutations are 
sparser, have higher noise resilience and fewer fixed 
point attractors as compared with networks outside of 
this wildtype component. 
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Neutrality [T is crucial for robustness and evolvability 
[2] of biological systems. It describes the fact that the 
mapping from genotypes to phenotypes is not invert- 
ible. A given phenotype can be encoded by more than 
one genotype. As Wagner [2 writes, ^'most problems 
the living have solved have an astronomical number of 
equivalent solutions, which can be thought of as existing 
in a vast neutral space". 

Computational studies of biopolymers revealed the 
existence of neutrality in the relation between sequence 
and spatial structure. RNA molecules and proteins are 
generated as a chain (sequence) of nucleic bases and 
amino acids respectively. The number of sequences fold- 
ing into one and the same functionally relevant spatial 
structure is found to be large. It is growing exponen- 
tially with the size of the molecule [H |4] . Together with 
an adjacency given by single mutations, the phenotypi- 
cally equivalent genotypes form the neutral network (or 
neutral graph). The properties of this graph, in partic- 
ular its connectivity, determine the robustness of the 
given genotype under mutations and its evolvability to- 
wards new phenotypes. 

Going from single molecules to the level of the whole 
organism, the phenotype is not given by the set of its 
molecule structures alone: The dynamics that arises as 
the result of activating and suppressing interactions be- 
tween molecules is crucial. This set of interactions is 
captured as a regulatory network [5] and gives rise to 
a temporal sequence of chemical concentration vectors 
that are responsible for the division of a single cell or the 
development of an embryo. Again, the mapping from 
genotypes (interaction networks) to phenotypes (tem- 
poral sequences) is not injective, i.e. several network 
structures are able to produce the regulatory dynamics 
of a given phenotype [H [7] . 

Here we apply the neutral graph concept to a dy- 
namical model [8 of cell cycle regulation in the organ- 
ism of the yeast species Saccharomyces cerevisiae (bud- 
ding yeast). In section[2]we introduce the model dynam- 
ics and the wiring of the wildtype network. The ensem- 
ble of functional networks that yield dynamics equiv- 
alent to the wildtype is analyzed in section [3) finding 
the neutral graph to be disconnected. In section [4] we 
focus on statistical properties of the subset of networks 
that are reachable from the wildtype. After a remark 
(section p| on the computation of network statistics, 
section |6] offers a discussion and open questions. 
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2 Cell cycle network and Boolean dynamics 

During the process of cell division, a eukaryotic cell 
grows and divides into two daughter cells. A cell cycle 
consists of four distinct and separate phases named Gi , 
5, G2 and M. In the Gi ("growth") phase, the cell com- 
mits itself for cell division under certain conditions. In 
particular, a necessary cell size must have been reached. 
A copy of the genetic information is produced in the S 
("synthesis") phase. The G2 ("gap") phase precedes the 
actual cell division in the M phase ("mitosis"). 

Here we are interested in the network of molecules 
(cyclins, inhibitors and degraders of cyclins and tran- 
scription factors) regulating this process. We consider 
the regulatory network of the monocellular eukaryotic 
organism Saccharomyces cerevisiae (budding yeast). Its 
genome comprises 13 million base pairs and 6275 genes, 
of which approximately 800 are involved in the cell cycle 
dynamics [9 . The dynamics is controlled by a core of 11 
key regulators with 34 directed interactions [8 , shown 
in Figure [2?TJ a), which we denote as the wildtype net- 
work. Interactions are captured by a matrix A. If node 
j has an activating effect on node the corresponding 
matrix element is a^j = +1, while inhibition is coded as 
aij = —1. In case of no direct influence from j to we 
have aij =0. Li et al. [8] model the regulatory dynam- 
ics with a Boolean approach [TOl [11] where each node i 
takes state values Si{t) G {0, 1} when being inactive / 
active at time t. In the time-discrete dynamics, nodes 
are updated synchronously, based on their weighted in- 
put sum hi{t) = aijSj{t). The state at the next time 
step is obtained by applying the threshold update rule 



1, hi{t)>0 
0, hi{t) < 
[Si{t), hi{t) = 



(2.1) 



From an initial condition 6'(1), representing the real 
starting state of the cell cycle, the dynamics produces 
the sequence of state vectors 5'(1), *S'(2), . . . , 6'(13), shown 
in Figure [2?TJb). The state 5'(13) = Gi is a fixed point 
of the dynamics. The system remains in this state until 
node Cln3 is externally activated. In the real system 
the external activation indicates that the cell size is 
sufficient for another division. 



3 Functional networks and the neutral graph 

Broadening our treatment of regulatory networks, we 
consider the set of all networks with interaction matri- 
ces over 11 nodes with entries a^j G { — 1,0,+1}. We 
call a network functional if it produces the state transi- 
tions of the cell cycle sequence in Figure [2?Tlb). Thus, 
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Fig. 3.1 Histograms of the number of interactions over func- 
tional networks. Positive (green curve), negative (red curve), 
curve) and total connections (black curve) of almost all func- 
tional networks exceed the corresponding counts in the wildtype 
network (vertical dashed lines). 



the wildtype network is functional. However, there are 
further functional networks. Out of the set of all 3^ ^ 
5.4 X 10^^ networks, approximately 5.11 x 10^^ are func- 
tional [12 . Figure [2?TJc) shows an example of a func- 
tional network different from the wildtype. 

Figure |3.1| shows the statistics for the number of 
interactions (arcs) present in functional networks. The 
wildtype network is sparse in comparison with the av- 
erage functional network. However, there are functional 
networks that are even sparser than the wildtype. These 
findings analogously hold when activating and inhibit- 
ing interactions are counted separately. Interestingly, 
functional networks have generally more suppressing 
than activating interactions, as is the case for the wild- 
type. 

A structure to refiect mutations on the set of func- 
tional networks is the neutral graph. Its nodes are the 
functional networks. Functional networks A and B are 
adjacent (connected by an edge) in the neutral graph if 
A is turned into 5 by a single mutation. According to 
our definition a mutation is a replacement of one entry 
in the interaction matrix. The Hamming distance be- 
tween two networks is the number of entries in which 
their interaction matrices differ. In order to avoid con- 
fusion with the networks of interaction we employ the 
term neutral graph as a synonym for the more com- 
monly used neutral network. An important property of 
a neutral graph is its connectedness. A mutational walk 
from network A to network 5 is a sequence of single 
point mutations that turns A into B without passing 
through non-functional networks. The neutral graph is 
connected if such a mutational walk exists for each pair 
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Fig. 2.1 (a) The Cell Cycle Network of the yeast wildtype has 11 nodes connected with activating (green) and inhibiting (red) 
interactions. Self-suppression is indicated by yellow loops, (b) A sequence of 13 states defines a cell cycle, as produced by the network 
in (a), (c) A different network (mutant) performs the same sequence of states. As the wildtype, this mutant has 34 interactions. 
However, 19 entries in the interaction matrix differ from the wildtype. 



of functional networks. We find that the neutral graph 
considered here is disconnected. One cannot pass from 
all functional networks to all others by sequences of 
mutations that preserve functionality. In fact, mutual 
reachability between functional networks is rare. The 
neutral graph falls into 4.7 x 10^ connected compo- 
nents with sizes distributed between 6.1 x 10^^ and 
^ 4.4 X 10^^, as shown in Figure 3.1 The component 



of the wildtype comprises around 5.66 x 10^^ functional 
networks. 
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Fig. 3.2 Cumulative size distribution of connected components 
of the neutral graph (falling curve). The component containing 
the wildtype has size 5.66 x 10^^ (vertical line). 



In this section we extend the analysis of the neutral 
graph. We focus on a comparison between functional 
networks in the wildtype component and all functional 
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Fig. 4.1 Comparison of statistics between all functional networks (green curves) and functional networks in the wildtype component 
(blue curves) of the neutral graph, (a) histogram of negative interactions, (b) histogram of positive interactions and (c) histogram of 
total number of interactions in functional networks, (d) histogram of Hamming distances (minimal number of mutations) from the 
wildtype. (e) Distribution of basin sizes of the Gi fixed point. The inset shows a zoom into the histogram for very large basin sizes. 
Histograms in panels (a)-(d) are exact. Histograms in (e) were obtained by uniform sampling of 10^ functional networks each from 
the whole neutral graph and from its wildtype component, respectively. 
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networks. Figure 4.1 ^a-c) shows how the number of (a) 
negative, (b) positive and (c) ah interactions is dis- 
tributed. Ah three plots reveal a significant statistical 
difference between networks in the wildtype component 
and the set of all functional networks. Networks in the 
wildtype component are sparse compared with the av- 
erage functional network. 

Geometric information of the neutral graph is pro- 
vided in Figure |4.1[ d) in terms of the Hamming dis- 
tance of functional networks from the wildtype. Func- 
tional networks in the wildtype component are closer 
to the wildtype than the average functional network is. 
Still the most remote networks in the wildtype compo- 
nent are found at distance 77 from the wildtype. Despite 
its moderate size, the wildtype component pervades a 
large part of the network space. 

Shifting attention from the structural to the dynam- 
ical properties of the functional networks, let us analyze 
the resilience of the dynamics against perturbations. As 
a measure of resilience we use the Gi basin size [8], i.e. 
the number of states from which the dynamics eventu- 
ally reaches the fixed point Gi. Clearly, the basin con- 
tains at least the 13 states in the cell cycle sequence. As 
shown by the distributions in Figure [XTJe), actual Gi 
basin sizes in functional networks contain many more 
states. Compared with all functional networks, basin 
sizes of networks in the wildtype component concen- 
trate at higher values. The most frequently observed 
basin size is 2047 for networks in the wildtype compo- 
nent, cf. the inset of Figure [XTj ^e). However, we have 
not found a functional network where the GI basin con- 
tained all 2048 states. Moreover, the distributions in 
the number of fixed points of functional networks show 
a striking difference between the wildtype component 
and the whole neutral graph. Figure [42| a) displays ge- 
ometric distributions in both cases. However, networks 
in the wildtype component have a significantly narrower 
distribution of fixed points. Interestingly, dynamic at- 
tractors (limit cycles) with more than one state show 
practically the same statistics in the wildtype compo- 
nent as in the whole neutral graph, cf. Figure 



4.2 



b). 



5 Computational aspects 



As noted by Lau et al. the set of network matrices 
performing a given state sequence has a simple com- 
binatorial structure. One can check independently for 
each node i if it takes the required state at each time 
step t. The states taken by node i only depend on the i- 
th row and not on the whole matrix. Thus, a functional 
network can be constructed by independently combined 
functional row vectors into a matrix. The set of func- 
tional row vectors for each node i is obtained by testing 



each of the 3^^ ^ 2x 10^ possible vectors over { — 1, 0, 1}. 
For observables, such as the number of interactions and 
Hamming distances that fall into sums over row vectors 
of the matrix, exact distributions are obtained by com- 
bining the distributions for all rows. In fact, each row 
of the matrix has its own neutral graph. The Cartesian 
product [13] of these is the neutral graph of the whole 
system. Sampling is used to obtain statistics of observ- 
ables that are not a function of single rows, such as the 
number of attr actors and basin sizes. 



6 Discussion & Outlook 

We have analyzed the neutral graph (also called neu- 
tral network) of discrete regulatory networks reproduc- 
ing the cell cycle sequence of budding yeast [8]. The 
neutral graph falls into many connected components. 
Networks in different components of the neutral graph 
are not accessible to each other through a sequence of 
mutations that retains cell cycle functionality. Our find- 
ing contrasts with the connected neutral graphs in the 
work by Ciliberti et al. in a similar type of discrete reg- 
ulatory networks [6[ |7]. There, function is defined as 
the eventual arrival at a predefined fixed point from a 
given initial condition. In the present study, the exact 
sequence of states leading to the fixed point is part of 
the required phenotype. We hypothesize that the frag- 
mentation of the neutral graph is caused by increasing 
functional constraints. 

Further analysis has revealed that functional net- 
works accessible from the empirical wildtype are struc- 
turally and dynamically distinct from other functional 
networks. Networks in the wildtype component are more 
sparsely wired and their dynamics is more resilient to 
perturbations, as compared to the average of all func- 
tional networks. 

Thus, networks in the wildtype component have prop- 
erties similar to the wildtype itself. This is remarkable 
since most networks in the wildtpye component are dis- 
tant from the wildtype, having only a few interactions 
in common. 

Future investigations could establish conditions for 
the connectedness of the neutral graph. To what extent 
is the fragmentation of the neutral graph caused by the 
strong discretization of interaction strengths? Allowing 
finer adaptations would lead to less fragmented neutral 
graphs. In the extreme (though chemically unrealistic) 
limit of continuously evolving interaction strengths, the 
set of all functional network matrices is convex and thus 
connected. 
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Fig. 4.2 The number of attractors of all functional networks (green curves) and the functional networks in the neutral graph 
component containing the wildtype (blue curves), (a) Distribution of the number of fixed points, (b) Distribution of the number of 
limit cycles (attractors of length at least 2). The wildtype itself has 7 fixed points (vertical dashed line) and no limit cycles. 
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